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Abstract. We present a Vlasov, i.e. a kinetic Eulerian simulation study of nonlinear 
collisionless ion-acoustic shocks and solitons excited by an intense laser interacting with 
an overdense plasma. The use of the Vlasov code avoids problems with low particle 
statistics and allows a validation of particle-in-cell results. A simple, original correction 
to the splitting method for the numerical integration of the Vlasov equation has been 
implemented in order to ensure the charge conservation in the relativistic regime. We 
show that the ion distribution is affected by the development of a turbulence driven 
by the relativistic “fast” electron bunches generated at the laser-plasma interaction 
surface. This leads to the onset of ion reflection at the shock front in an initially cold 
plasma where only soliton solutions without ion reflection are expected to propagate. 
We give a simple analytical model to describe the onset of the turbulence as a nonlinear 
coupling of the ion density with the fast electron currents, taking the pulsed nature of 
the relativistic electron bunches into account. 


PACS numbers: 


Submitted to: Plasma Phys. Control. Fusion 

1 . Introduction 

The production of high energy protons by laser-driven acceleration is currently a subject 
of active research. Several experiments have demonstrated the generation of multi-MeV 
proton beams in a wide range of laser and target parameters Different mechanisms 

have been studied in order to optimize and control the main characteristics of the proton 
beam for the foreseen applications. 

Among the various ion acceleration mechanisms, here we focus on collisionless shock 
acceleration (CSA). The simplest CSA scenario may be briefly described as follows. 
Collisionless shock waves can be excited at the laser-plasma interaction surface due to 
the combination of radiation pressure-driven pushing and steepening of density prohle 
and of electron heating up to high temperatures. In the electrostatic regime (to which 
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we restrict ourselves in the present paper) the shock waves are characterized by an 
electrostatic potential jump at the shock front, moving with velocity Vg. Some ions can 
be accelerated from rest by reflection from the shock front as a moving wall, with a 
velocity gain of 2Vs, provided that the ions have a sufficient velocity along the direction 
of the shock propagation so that they cannot overcome the potential barrier. For a 
plasma with electron temperature Tg, the shock velocity is 14 = Mcg where M > 1 is 
the Mach number and Cg = (ZTf,/mi)^'^^ is the ion acoustic velocity with mi = Amp is 
the ion mass, rup the proton mass, Z and A the ion charge and mass numbers. At laser 
irradiances > 10^® W cm“^/im (with I and A the laser intensity and wavelength, 
respectively) “fast” high-energy electrons are generated with effective values of of 
several MeVs. Thus, reflection from the shock front may generate multi-MeV ions with 
typical energy Si = mi{2Vg)‘^/2 = 2M^ZTe, and a monoenergetic spectrum as far as the 
shock velocity remains constant. 

The CSA mechanism has been explored previously via particle-in-cell (PIC) 
simulations [Sflll] and in recent laboratory experiments using linearly polarized, CO 2 


laser (A = 10 /rm) pulses 12-15 and gas jet targets which were slightly overdense 
for the laser pulse, i.e. having an electron density Ug > ric = 47rmeC^/(e^A^), the 


cut-off (or “critical”) density. In the experiment by Haberberger et ah 12 


a narrow 


monoenergetic peak corresponding to ~ 20 MeV protons was observed using a linearly 
polarized train of pulses with duration of a few picoseconds. While the possibility to 
accelerate monoenergetic protons is appealing, the proton flux observed in Ref. 12 is 
relatively low (~ 10^ MeV~^sr“^) so that the population of accelerated protons would be 
hardly resolved in a PIC simulation, unless the number of particles per cell is drastically 
increased up to values that would be computationally challenging for multi-dimensional 
simulations. Moreover, PIC simulations of CSA have shown the crucial role of the 
background ion temperature on shock formation and particle reflection [^, so that the 
high resolution of the low-density tail and of the non-thermal component in the ion 
distribution function is important. This issue suggests that PIC simulations should be 
tested when possible against Eulerian “Vlasov” simulations which are free from effects of 
low particle statistics and fluctuations, at the cost of higher computational requirements. 

In this paper, we studied CSA using a Vlasov code and comparing the results 
to those from PIC simulations. The numerical approach is described in section The 
Vlasov code includes a simple method to adapt the splitting algorithm to the relativistic 
case without violating mass conservation, as described in the Appendix. Our results 
include the observation of acceleration at the shock front also in the case with no thermal 
spread for the initial ion distribution function (Tj ~ 0) where the theory (resumed in 
section predicts only the generation of a solitary wave, i.e. an ion-acoustic soliton, 
without reflection from the shock front. This observation is related to the development 
of a sort of ion density turbulence in the upstream region ahead of the shock front, 
locally broadening the ion velocity space and allowing for ion reflection. In section 
we give a simple analytical model in which, on the basis of the observations, we relate 
the ion turbulence to the pulsed nature of the fast electron bunches generated by the 
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2. Background theory 


The formation of a collisionless electrostatic shock solntion is not possible in principle 
within the flnid description, where only soliton solntions appear as the result of a 
balance between nonlinearity and dispersion. A shock wave solution may exist in 
presence of kinetic effects acting as an effective dissipation, breaking the symmetry 
of the soliton 16 . In this framework, ion acceleration due to reflection from the shock 


front may be a mechanism inherent to the collisionless shock formation rather than a 
consequence of it. 

A theory of collisionless nonlinear ion-acoustic shocks and solitons is resumed in 
Refs. 16,17 where the electrons are assumed to be non-relativistic and in Boltzmann 


equilibrium. In the shock front rest frame, moving at velocity Vs with respect to the 
laboratory frame, the reflection occurs for all the ions whose kinetic energy does not 
exceed the height of electrostatic potential barrier (fm- All the ions in the high energy tail 
of the distribution function with a velocity component Vi along the shock propagation 
direction in the laboratory frame such that 


Vi> Vs- 



( 1 ) 


will get reflected from the shock. Thus, a spread in the initial velocity distribution 
function is required, otherwise no reflection occurs and only soliton solutions are possible. 
The number of reflected ions should be however small enough to avoid loading of the 
shock, e.g. excessive energy transfer to the accelerated particles [^. This makes clear 
the importance of the initial ion distribution in order to predict the development and 
lifetime of a soliton or a shock wave and the fraction of accelerated particles. 

Heuristically, the formation of high-speed shocks or solitons in intense laser 
interaction with overdense plasmas involves two requirements: the electrons must be 
heated to high temperatures to allow shock/soliton propagation, and a strong initial 
density perturbation must be driven. For what concerns the hrst requirement, a linearly 
polarized pulse generates a fast electron population with a typical kinetic energy 

Tp = ruec^ ((l + ao/2)^/^ - l)) , (2) 

where Oq = ifll = 0.85(/A^/10^® W cm is the peak dimensionless 

amplitude of the laser held. If we roughly assume the fast electrons to have a 
Boltzmann distribution with temperature Tj ~ Tp, then the fast electrons can sustain 
the propagation of nonlinear waves characterised by the ion-acoustic velocity Cs = 
[ZTj. Recirculation of the fast electrons through the target is also necessary 
to produce a uniform temperature, which may introduce a dependence upon the pulse 
duration and target thickness. 

The driving perturbation can be provided by the radiation pressure of the laser 
pulse Pl = {1 + R)I/c (with R < 1 the rehection coefficient) pushing the plasma 
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surface as a piston with an average velocity Uhb, commonly known as the hole boring 
(HB) velocity. The expression of Uhb can be obtained either via a momentum balance 
equation 18 or via a dynamic model as (for Uhb c) 


Uhr — 


Pr. 


2miP 


1/2 


1 -\- R Z rUe ric 
2 A mp fit 


1/2 


Uo 

71 ' 


( 3 ) 


The dynamic model shows that a fraction of ions is accelerated up to a velocity 2 uhb 
( following wave breaking of the density profile (^). The generation of such ions is also 
needed for momentum conservation 18,20 . Density perturbations may detach from the 
surface and propagate with supersonic speed (M > 1) into the plasma when Uhb ^ Cg. 

A collisionless shock may also be generated due to the nonlinear evolution 
of instabilities, which in the multi-dimensional case include electromagnetic 
counterstreaming (or “Weibel”) instabilities driven by fast electrons 21 . In Ref. 12 
the authors advocate some driving mechanism different by radiation pressure in order 
to account for a shock velocity Vg 3> Uhb- In the present paper we restrict to one¬ 
dimensional (ID) geometry and focus on radiation pressure-driven shocks. 

Notice that in a cold plasma, i.e. in the absence of fast electrons, a “true” shock 
or soliton wave may not form and detach from the surface. In this case, only the 
above mentioned ion population with velocity ~ 2 uhb is observed; such ions move into 
the plasma ballistically, as charge-neutralized bunches 19 . This is the so-called HB 


acceleration regime. Although one may consider such ions as “reflected” from the laser- 
driven piston, so that the process sounds quite similar to CSA, there are basic differences 
between CSA and HB acceleration, which is favored if circular polarization (and normal 
incidence) of the laser pulse is used instead of linear polarization since in the former 
case the fast electron generation is strongly reduced 19 . While the shocks driven in 


CSA accelerate ions along the propagation in the plasma, HB occurs only at the plasma 
surface and during the laser pulse action. In addition, the number of ions accelerated 
by HB is larger (the density of ion bunches may be close to the background value). This 
might explain why the flux of detected protons reported in Ref. 12 is some five orders 
of magnitude lower than that observed with a similar experimental set-up, but using 
circular polarization 22 . In this case the proton spectrum shows a broader peak at the 


lower energy ~ 1 MeV which is fairly consistent with HB theory. In the case of linear 
polarization, HB and CSA may both occur. 


3. Numerical method and simulation set-up 

We developed a Vlasov code that provides a solution of the ID Vlasov-Maxwell system 
of equations within a completely Eulerian approach. Once the distribution function of 
each species is discretized on a fixed grid in the real space and in the momentum space, 
the code performs a direct integration of the Vlasov equation. The Vlasov equation is 
coupled with Maxwell’s equations for self-consistent electromagnetic fields, generated 
by the charge densities and currents of the particle species. These latter quantities 
are defined as the momenta of the distribution function and are calculated over all the 
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Eulerian grid cells. In the absence of a magnetic field component along the propagation 
direction = 0), the Vlasov eqnation can be rednced to a IDIP geometry exploiting 
the conservation of transverse canonical momentnm fl^ = p± + (lA±, where A± denotes 
the transverse component of the vector potential. Thus for a particle species with mass 
m and charge q, the evolution of the reduced one-dimensional distribution function 
/ = f{x,px,t) is described by 


dt my dx 


-V q \ Ex 


q \ df 

dpx 


2m7 


dx 


= 0 


(4) 


For the numerical integration of eq.Q, the Time Splitting Scheme 23 and the Positive 
and Flux Conservative Method 24 have been employed along with an original method 
providing the exact mass conservation in the relativistic regime (see Appendix A). 


Violation of mass conservation was a known drawback of the use of splitting schemes in 

In order to ensure mass conservation quite complex methods 
Our method is much simpler than others, but 


the relativistic regime 25 


have been previously proposed 25,26 
apparently effective. 

In the next section, we report the results of a simulation with an electron-proton 
plasma of length L = 10 A having a steep density profile with a linear rising and falling 
ramp of length 0.2 A, where A is the laser wavelength. The plasma plateau density 
is set to Hi = He = 2.0 ric- The laser pulse has linear polarization, peak amplitude 
Oo = 2.0 and duration r = 60 T, with T the laser period. The temporal prohle has 
a sin^-like rising and falling ramp of one period length, and 58 T of constant plateau. 
The temporal and spatial resolution is set to Ax = cAt = A/10^. The size of the 
simulation box is = 16 A. In the momentum grid the resolution is Ap^^e = 0-05 rUeC 
and Apx^i = 0.045 mgC, respectively for electrons and ions. The distribution function, 
for each species, has a Gaussian shape in momentum space with initial temperatures of 
Te = 5 KeV and T* = 1 eV. A lower electron temperature would require a much greater 
computational effort, and as soon as the electrons are heated up by the laser-plasma 
interaction, the spread of the distribution function becomes much higher than the initial 
distribution. A slightly higher resolution is chosen for the ion momentum with respect 
to the electrons, in order to resolve the initial spread of the ion distribution function 
that is initialized with a temperature much lower than the electron distribution. Open 
boundaries are used for the fields and reflecting boundaries for the distribution function. 
We verified that the density of particle reflected at the boundaries is small enough to 
have a negligible role on the physical results. 


4. Simulation results 

4-1. Vlasov simulation of shock generation 

Figure shows snapshots of the ion density Uj, the electric field E^, and the ion phase 
space at different times. The incoming laser pulse propagates from left to right. At the 
first time shown in figure = 17T), the sharp density peak generated by the radiation 
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pressure push is apparent. The density peak moves with a velocity of Vp ~ 0.026c, as 
given by the fit in Figj^ which is slightly higher than the value of Uhb — 0.023 (since 
R ~ 0.9 in the simulation) predicted by Eq.(|^. Our explanation is that, in this regime 
(low density and linear polarization) a significant number of electrons are dragged on 
the vacuum side, as it is apparent from the phase space plot of the electrons in figure 
This causes some ions to be accelerated in the backward direction, i.e. towards the 
incoming laser pulse, as observed in the ion phase space plot (figure [^. Since the total 
momentum flux from the laser to the plasma is always given by Pl and it is in the 
forward direction, the ions accelerated forward must get an extra boost in addition to 
that leading to ([^. In other words, fast electrons produce a sort of “ablation pressure” 
which adds up to the laser radiation pressure. 

The phase space plot at t = 17T in figure shows that the ion density profile is 
undergoing wavebreaking, which at t = 27T has produced to the typical “X”-type phase 
space structure 27-29 with further acceleration of some ions which cross the breaking 
point 19 . Correspondingly, the density peak appears to split up in two peaks, the 


second moving at constant velocity of Vs = 0.039c as also shown by the fit in figure [T] 
At later times, an ambipolar electric field structure is evident around the density peak, 
while a “pencil” of ions reflected at velocity ~ 2Vs appears in the phase space. We 
thus identify the fast density peak as the front of a radiation pressure-driven shock. 
A fit on the electron spectrum gives a fast electron temperature of ~ O.SdrUeC^ (to be 
compared with Sp = 0.73meC^ from Eq.([^, which yields Cg — 0.021 and a Mach number 
M = Vs/cs ^ 1.9. 

The high energy or “fast” electrons produced by the interaction of the high-intensity 
laser with the front surface of the target propagate through the bulk as small-duration 
bunches (see the electron phase space plot in figure]^. When the electrons reach the 
rear side of the plasma they produce a negative charged sheath. Because of the electric 
field generated by the expanding sheath at rear side, most of the electrons cannot escape 
and carry on recirculating across the target. The recirculation leads to a spread of fast 
electrons all across the target. As shown in figure the first recirculating electrons 
reach back the front side around t = 27T, later than the formation of the shock front. 
The latter keeps a constant velocity also at later times, as shown in figure Thus, 
we conclude that the shock is driven by the radiation pressure action and the related 
wavebreaking, as the onset of fast electron recirculation does not affect the velocity of 
the shock. 


4 . 2 . Ion acceleration and turbulence 

The density and field profile at t = 82T in Figj^ show two prominent density peaks 
associated to ambipolar field structures. A third, much weaker structure is also visible. 
The propagation velocity of the second density peak is slower than the first one. This 
latter observation may suggest that the structures may be better described as two 
solitary waves, i.e. ion-acoustic solitons, rather than a coherent shock wave. Indeed, as 
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Figure 1. Ion density and longitudinal electric field (in units of Eq = meCtole) (top 
row) and the ion distribution function in the phase space x — Px (bottom row) at time 
t = 17, 27, 82 in units of the laser pulse period T = A/c = 2'kIuj. The inset shows the 
position vs. time of the rightmost density peak, which is fitted by two linear functions. 




Figure 2. Electron distribution function in the phase space x — Px at times t = IIT 
and t = 27T for the Vlasov simulation, and at times t = 27T for the PIC simulation. 


predicted by theory 16,17 , in the case of cold ion population we actually do not expect 
a shock wave but a soliton solution to propagate. For this solution the propagation 
velocity increases with the amplitude of the electrostatic potential associated with the 
density peak. However, the ion phase space plot for t = 82T in hgure clearly shows the 
presence of ions reflected from the hrst propagating structure, which is a characteristic 
feature of shock solutions 16 10 (while the “curly” feature behind the front corresponds 
to ions trapped between the two structures). 

Figure l^shows the ion energy spectrum at time t = 82T, which has peaks at energy 
S ~ 4.7meC^ and E ~ 5.7meC^ — 2.9MeV with energy spread ~ 2%. The peak at higher 
energy corresponds to the plateau with momentum ~ 145 nieC in the phase space 
(£gure[^. This is very close to twice the propagation momentum of the hrst density 
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Figure 3. Ion energy spectrum at t = 82 T in the x-range 6.3A ~ lOA for the Vlasov 
(left) and the PIC (right) simulations, respectively. 


peak P ~ 72.1 mgC, as expected for the reflection from a shock front. 

The second energy peak aronnd T ~ 4.7 rUeC^ corresponds to a gronp of ions which 
are ahead of the first density peak at a: ~ 8A and have thns been generated previonsly, in 
an early stage of radiation pressnre pushing and shock formation and before the stable 
shock propagation regime. Overall, the ions contained in the two spectral peaks contain 
a fraction < 10“^ of the laser pulse energy. 

In order to estimate the number of reflected particles according to eq.([^, we 
measure the value 0^ of the electrostatic potential jump around the shock front, i.e. 
the fastest density peak. For the observed value ecfm — 1.3 the reflection condition 
([^ is fulhlled by the ions with velocity Vi ~ 70 Vth, where Vth is the thermal velocity 
corresponding to the initial ion temperature T* = 1 eV. In such region of the phase 
space the tails of the distribution function are completely negligible (~ e“^°^), thus we 
should not expect to observe reflected particles, in contrast to what is obtained from 
the simulation. 

The presence of accelerated particles is justihed by the growth of a perturbation in 
the ion density in the upstream region, where the shock has not yet propagated, as shown 
in figure]^ a). In correspondence with the fluctuations of the ion density, an oscillation 
of the ion distribution function is observed in the phase space, see figure |^b). The 
oscillations of the ion velocity produce a spread of the distribution function and allow 
a fraction of ions to exceed the threshold value of Eq.([^, so that they get reflected by 
the shock front. The turbulence in front of the first ion density peak leads to a variable 
quantity of reflected ions with time. 

In order to characterize the ion density perturbation we notice that it starts growing 
at the rear side of the plasma when the fast electrons have passed through this region, 
after being pushed back inside the target by the electrostatic held of the sheath, and that 
it has a regular periodicity Aj ~ A/4, as shown in £gure|^a). These observations will 
guide us to sketch a simple model for the generation of the ion density perturbation, 
in which the temporal structure of the fast electron bunches plays an essential role 
(section 1^. 
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Figure 4. Left frame: growing oscillations of the ion density showing a spatial period 
equal to A/4, with A the laser wavelength. Right frame: ion distribution function in 
the phase space x — Px and ion density (red line) at t = 67 T 


4.3. Comparison with PIC simulations 

In this paragraph we show the comparison between the results obtained with the PIC 
approach and the ones obtained with the Vlasov code. With respect to the Vlasov 
approach, the well known drawbacks of PIC simulation are the high statistical noise 
level and the limited resolution of low density regions of the phase space. In order to 
test the convergence of the PIC simulations, runs were performed for three different 
numbers of particles per cell per species Np = 10^, 10^, 10^ and also for two different 
values of the spatio-temporal resolution Ax = cAt = A/10^ and Ax = A/10^. 

In £gure[^we compare the Vlasov results on the development of ion turbulence with 
the PIC ones obtained with Np = 10^ and the two different choices of Ax. Despite the 
relatively high resolution, the noise level in PIC simulations makes it difficult to clearly 
characterize the ion perturbation growth and to accurately measure the periodicity of 
the oscillations, whereas at later times the noise leads to a turbulence in the upstream 
region with greater amplitude than in the Vlasov case. The difference can be attributed 
to the higher level of initial noise which acts as a seed for the nonlinear three-wave 
process generating the turbulence (section]^. It may also be noticed that in the Vlasov 
code the density perturbations are smoothed to some extent by numerical diffusion. 

The difference in the amplitude of the ion turbulence affects the ion reffection from 
the shock front since the broadening of the distribution function in the upstream region 
is strictly dependent on the characteristics of the turbulence. The different fraction of 
particles, in the two numerical approaches, lying in a region of the phase space where 
condition ([^ is fulhlled, may vary the transfer of kinetic energy from the shock front to 
the accelerated particles. As already discussed in Ref. [^, a high number of accelerated 
particles may slow down the shock and at later times may broaden the energy spectrum 
peak toward lower energies. In £gure|^the ion phase space from the PIC simulation at 
t = 82T is superimposed to the ion phase space obtained with the Vlasov code. The 
reflection in the PIC simulation is not as steady as in the Vlasov case and the shock 
velocity is slightly higher. Differences are also apparent in the energy spectrum of the 
accelerated ions (figure [^, with a third spectral peak appearing at ~ 7 MeV. Thus, the 
comparison show that the PIC simulation predicts higher energy and efficiency. Finally, 
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Figure 5. Ion density oscillations aXt = 27T (left frame) and t = 77T (right frame). 
PICl corresponds to the simulation with Ax = A/10^, PIC2 with Ax = A/10^. 



Figure 6. Ion distribution function at t = 82r for the Vlasov and PIC (purple dots) 
simulation with Ax = A/10^ (left frame). 


a comparison is also shown for the electron phase space in hgure For this latter, the 
very good agreement is a test of the accuracy of the relativistic splitting scheme for the 
Vlasov code, including a correction for mass conservation (see 


Appendix A). 


5. A model for ion turbulence 


In this section we present a simple model which accounts for the generation of ion density 
oscillations with a spatial periodicity of ~ A/4, where A is the wavelength of the laser 
pulse. 

In the case of a linearly polarized intense laser pulse incident on an overdense 
plasma, electrons are accelerated by the magnetic component of the Lorentz force 
(J X B) and are pushed into the plasma as bunches at a 2u rate, being oj the laser 
frequency. Here we focus on the most energetic electrons propagating at velocities close 
to c for high enough intensities. These fully relativistic electron bunches maintain their 
coherence during the propagation. An experimental proof of such coherence comes from 
the spectrum of the optical transition radiation emitted when the electrons reach the 
rear boundary of the target 30 . The lower energy part of the distribution function 


corresponding to velocities < c will lose temporal coherence. In our model we thus 
assume the coherent, fully relativistic electrons as a population separated from the 
background electrons which have been heated up forming a Boltzmann distribution 
with a typical temperature Tg, assumed to have a non-relativistic value. 
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Figure 7. Schematic representation of the relativistic electrons bunches propagating 
in the positive direction (blue) and of the bunches reflected by the electric field in the 
expanding sheath (red). 


Most of the relativistic electrons reflux back inside the plasma because of the 
electrostatic field that acts as a reflecting wall. Thus, after a time interval which depends 
on the target size there are two counterstreaming populations of relativistic electrons, 
both bunched with spatial periodicity Xp ~ 27ic/{2u) = A/2 as it can be noticed in 
hgure The situation is sketched in hgure 

In our model we consider the nonlinear beat of the counterpropagating relativistic 
bunches, which we consider as two pump waves for the system. We write the electric 
held perturbation associated to the bunches as 

EP{x, t) = ^e^k+x-iu^+t ^ ^ ^ ^ ^5^ 


with u:± ~ 2a; and k± k. 27r/(Ap) = dvr/A. The resulting nonlinear force on the electrons 
at the beat frequency may be calculated as follows. Starting with the equation of motion 
in the held ([^, 

= -eEP{xe{t),t) , ( 6 ) 

we solve for Xe(t) by a standard iterative method. Writing Xe(t) = Xoe(t) + Xie(t), we 
have 


d^x 

- -eE'P{xe{V),t) , 


- -eXoed:^EP{Xe{0),t) 


(7) 


Thus, in the equation for Xie(t) there is a term at the beating frequency fl = 

and wavevector K = k+ + k_. The response at such frequency can be described by the 

nonlinear force 


/n 


f- { 'tE.e: 

Arue Vi^+ 


\e_EI\ + c.c. . 

CJl / 


( 8 ) 


Notice that /nl = 0 for a symmetric situation with Ej^ = E^, a;+ = U- and = /c_. 
In our case this would correspond to elastic rehection of the relativistic bunches at the 
rear plasma boundary. The symmetry is broken because of the loss of energy to the 
background electrons from the bunches while propagating in the plasma, and at the 
plasma boundary via the self-generated sheath held as the bouncing electrons enter 
the vacuum side. In particular, since the background electrons have a temperature 
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Tg, the plasma boundary will be in motion with a typical velocity of the order of 
Cg = the ion-acoustic velocity. The reflection of the relativistic electron 

bunches from the wall moving at velocity u ~ Cs leads to a shift in the reflected frequency 
that gives different amplitudes between the forward and backward perturbations and 
a non-zero value of = a;+ — a;_ ~ 2u{u/c) (since u c). In particular, ~ Kcg, 
so that the nonlinear perturbation may couple efficiently to ion-acoustic waves. Such 
perturbation has wavevector = fc_|_ -|- /c_ ~ 2k± = 27r/(A/4) which agrees with the 
periodicity observed in the simulations. 

The coupling of the nonlinear perturbation at the beating frequency to ion- 
acoustic waves may be described by the system of linearized fluid equations 


dtUe + nodxVe = 0 , 

o 6 Tg dxVle 

OtVe = - E - 

rUg mg no 

dxE = 4:7re{ne - Zui) , 

d„H + = 0 , 

a,vi = —E. 

rui 


+ /n 


(9) 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 


For further simplihcation we assume that on this low frequency scale the electrons can 
be considered in a mechanical quasi-equilibrium condition, Ug ~ 0, and quasi-neutrality 
may be assumed, Ug ~ Zni. In this way we obtain for the ion equation of motion 


O ry 2 

dtVi ^ -ZCg 


no 


Z 

S /nL 5 

nrti 


so that, using the equation of continuity to eliminate Uj we eventually obtain 


(5* 

with the solution 


Csdx)Vi = -^t/NL , 

rUi 


ikl 


Vi = 


+ c.c. 


(14) 


(15) 


(16) 


mi (fi2 _ 

where /nl is the complex amplitude of /^l [Eq.([^]. Eq.(16) shows a resonant behavior 
when the fast electron-driven perturbation frequency is close to that of the ion acoustic 
waves. 

It may be argued that the above described effect is strongly enhanced by the ID 
geometry as the reflected electrons are in the same directions as the incoming ones. 
In a more realistic multi-dimensional geometry the electrons would reflux with some 
angular spread and the modeling would be more complicated. Nevertheless, our picture 
shows that in principle the bunched nature of fast electrons, commonly neglected in 
the modeling of their transport through the plasma, may play an important role in the 
development of nonlinear effects and instabilities. 
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Generation of laser-driven collisionless shocks and related ion acceleration has been 
studied by means of numerical Eulerian simulations, based on a corrected splitting 
scheme for the relativistic Vlasov equation. The high resolution of low-density phase 
space regions and the absence of numerical noise typical of Vlasov simulations allowed 
a clear characterization of the shock formation process and the observation of a new 
mechanism of fast electrons-driven ion turbulence. Hence, the Vlasov simulation 
approach appears to be a useful tool for the study of collisionless shock acceleration 


physics and it is being used already by other groups 31 


In the regime we investigated, the shocks are driven by the radiation pressure 
action and related wavebreaking in the ion density prohle. Ion reflection from the 
shock front occur also for initially cold ions because of the velocity spread generated by 
the development of ion turbulence in the upstream region. Particle-in-cell simulations 
show a higher amplitude for the turbulent oscillations which leads to differences in the 
spectrum of accelerated ions with respect to Vlasov simulations. A simple model in 
which the ion turbulence is driven by beats induced by pulsed, coherent, relativistic 
electron bunches has been introduced. 

Appendix A. Mass conservation in the splitting scheme 


In our code the so called Time Splitting Scheme 23 is exploited in order to perform 


the numerical integration of the Vlasov equation. The scheme has been widely used and 
benchmarked in the electrostatic non relativistic case; in order to calculate the evolution 
of the distribution function it treats separately in the Vlasov equation the convective 
terms in the x direction and the one along the momentum axis, obtaining a scheme 
accurate up to the second order. 

In the fully relativistic electromagnetic case, it has been demonstrated that the 
scheme does not conserve the particle density [^. Indeed in our case the Vlasov 
equation can be splitted as 


d,h + 4 1 —/i) = -/i4 (— 


dtf2 + dp^ {FL,xf2) = -f2dpS- 




d I A, 


(A.l) 

(A. 2 ) 


2 m 7 dx 

where is Lorentz force along the a:—axis and 7 = (1 -|- -|- q^A\{x)]/m?‘Y^‘^. 

The presence of the relativistic 7 factor leads to the appearance of extra terms 
on the r.h.s., which are vanishing in the electrostatic non-relativistic case for which 
7=1 25 . Applying the method presented in Ref. to calculate the evolution of the 


distribution function, a cumulative systematic error would be introduced at each time 
step, resulting in a poor conservation of the charge density for each species. Without 
applying any correction, an unphysical loss of charge is found in the region where the 
electromagnetic effects are important (i.e. in correspondence of the laser pulse, where 
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Figure Al. Comparison of the longitudinal electric field at t = 75 T for a Vlasov 
simulation carried out with (blue line) and without (red line) the correction and a PIC 
(green line) simulation. 


7^ 0). On the other hand, in the region where the transverse vector potential vanishes 
the system is automatically reduced to the electrostatic case. 

In order to ensure the total particle mass conservation the quantities on the r.h.s. 
of Eq.(A.2) are considered as local source terms. Therefore these terms are integrated 
at each time step and the corresponding density is reintroduced in the distribution 
function reproducing in the momentum space the distribution of the remaining particles. 
As an example, in Fig jAl| we report the comparison between a simulation carried out 
without the correction and a corrected one, running with the same parameters. We also 
show a PIC simulation that provides a test of the correctness in the relativistic and 
electromagnetic regime of the correction developed in our Vlasov code. 

Fig,^ shows the longitudinal electric held associated with the ion density peak 
propagating in the overdense plasma. If the correction is not applied, the loss in the 
electron density causes the plasma to be non exactly neutral and an electric held arises 
at the right boundary. Reintroducing the lost electron mass, the held vanishes at 
the end of the box. 
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